! shock_daemon module
! Purpose:
!   1. detect collision between parcels
!   2. mix the collided parcels
! Author:
!   Name: DONG Li
!   Email: dongli@lasg.iap.ac.cn

module shock_daemon

    use space_manager
    use time_manager, only: dt
    use tree_struct, only: forest, num_tree

    integer, parameter :: max_num_mix = 100
    integer, parameter :: max_num_parcel = 100

    type mix
        integer :: n
        integer :: i(max_num_parcel)
        real :: mix_time
        type(mix), pointer :: prev
        type(mix), pointer :: next
    end type mix

    integer :: mix_counter
    type(mix) :: mix_record(max_num_mix)

    real :: min_distance = 1.0e-4

contains

    subroutine shock_daemon_detect_mix(n, x, u)
        integer, intent(in) :: n
        real, intent(in) :: x(num_dim,n)
        real, intent(in) :: u(n,num_dim)

        integer :: i, j, k, i1, i2
        real :: x1(num_dim), x2(num_dim)
        real :: u1(num_dim), u2(num_dim)
        real :: dx(num_dim), du(num_dim)
        real :: temp1, temp2, temp3, temp4
        real :: t1, t2

        mix_counter = 0

        do k = 1, num_tree(3)
        do j = 1, num_tree(2)
        do i = 1, num_tree(1)
            do i1 = 1, forest(1)%tree(i,j,k)%np
            do i2 = i1+1, forest(1)%tree(i,j,k)%np
                x1 = x(:,forest(1)%tree(i,j,k)%pi(i1))
                x2 = x(:,forest(1)%tree(i,j,k)%pi(i2))
                dx = x2-x1
                u1 = u(forest(1)%tree(i,j,k)%pi(i1),:)
                u2 = u(forest(1)%tree(i,j,k)%pi(i2),:)
                du = u2-u1
                temp1 = sum(du**2)
                if(temp1 == 0) cycle
                temp2 = sum(du*dx)
                temp3 = sum(dx**2)-min_distance**2
                temp4 = temp2**2-4*temp1*temp3
                if(temp4 < 0) cycle
                temp4 = sqrt(temp4)
                t1 = (-temp2+temp4)/(2*temp1)
                t2 = (-temp2-temp4)/(2*temp1)
                temp1 = min(t1, t2)
                temp2 = max(t1, t2)
                if(temp1 < 0 .and. temp2 > dt) then
                    cycle
                else if(temp1 < 0 .and. temp2 > 0 .and. temp2 < dt) then
                    mix_counter = mix_counter+1
                    temp3 = temp2
                else if(temp1 > dt) then
                    cycle
                else if(temp1 > 0 .and. temp1 < dt) then
                    mix_counter = mix_counter+1
                    temp3 = temp1
                end if
                if(mix_counter > max_num_mix) then
                    ! handle error
                    write(*, *) "mix_counter exceeds max_num_mix"
                    stop
                end if
                mix_record(mix_counter)%mix_time = temp3
                mix_record(mix_counter)%n = 2
                mix_record(mix_counter)%i(1) = i1
                mix_record(mix_counter)%i(2) = i2
            end do
            end do
        end do
        end do
        end do

        return
    end subroutine shock_daemon_detect_mix

end module shock_daemon
